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1.  INTRODUCTION 

This  is  the  Final  Report  for  the  grant  “Evaluation  of  Acoustic  Propagation  Paths 
into  the  Human  Head,”  Grant  Number:  FA9550-06-1-0128,  for  the  period  between 
December  1,  2005  and  November  30,  2008.  The  Air  Force  Program  Manager  is  Dr. 
Willard  Larkin,  and  the  Principal  Investigator  at  University  of  Illinois  at  Urbana- 
Champaign  is  Professor  William  D.  O’Brien,  Jr. 

The  overall  objective  of  the  program  was  the  further  development  and  validation 
of  an  acoustic  wave  propagation  model  using  well  understood  and  documented 
computational  techniques  that  track  an  air-borne  incident  acoustic  wave  propagated 
around,  into  and  in  the  human  head.  The  validation  testing  of  the  computational  model 
involved  human  subject  testing.  The  theoretical  model  provided  a  known  fundamental 
basis,  and  the  experiment  provided  the  evidence  to  evaluate  the  model.  The  goal  to 
further  development  and  validation  of  an  acoustic  wave  propagation  model  has  been 
achieved. 

In  summary,  the  highly  interdisciplinary  program  has  addressed  successfully  the 
challenging  technical  goals  set  forth  to  develop  and  validate  an  acoustic  wave 
propagation  model  using  well  understood  and  documented  computational  techniques  that 
track  an  air-borne  incident  acoustic  wave  propagated  around,  into  and  in  the  human  head. 
A  finite  element  model  was  formulated  for  simulating  the  acoustic  wave  propagation  into 
concentric  spheres.  This  model  was  verified  against  the  closed  form  solution  of  the  plane 
wave  propagation  into  concentric  fluid  spheres.  A  ray  tracing  method  was  developed  and 
verified  to  visualize  the  acoustic  wave  propagation  pathways  from  a  scalar  pressure  field. 
Finally,  the  occlusion  effect  provides  a  relatively  robust  phenomenon  to  compare  the 
computational  findings  with  the  behavioral  and  functional  findings;  the  occlusion  effect 
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contour  (or  difference  score)  between  the  open  and  occluded  ear  show  a  similar  trend 
between  the  computational  findings  and  the  behavioral  and  functional  findings.  In 
conclusion,  the  computational  approach,  utilizing  fluid  elements,  provides  very  good 
agreement  with  human  subject  findings. 

Acknowledgements:  Professor  Margaret  Wismcr  (Department  of  Electrical 
Engineering,  Bucknell  University  and  Department  of  Electrical  and  Computer 
Engineering,  University  of  Illinois  at  Urbana-Champaign)  successfully  developed  the 
computational  model,  and  closely  interacted  with  ECE  graduate  students  Jared  A. 
McNew  and  Alessandro  Beilina.  Professors  Charissa  R.  Lansing  and  Ron  D.  Chambers 
(Department  of  Speech  and  Hearing  Science,  University  of  Illinois  at  Urbana- 
Champaign),  in  association  with  SHS  graduate  students  Sarah  Melamed,  Lynn  M.  Brault 
and  Woojae  Han,  successfully  guided  the  human  subjects  testing  that  yielded  new  and 
interesting  results  relative  to  this  grant’s  goals. 


2.  ACOUSTIC  WAVE  PROPAGATION  MODEL 


This  section  presents  the  numerical  simulation  of  evaluating  pathways  of  bone 
conducted  sound.  Two  separate  algorithms  have  been  developed  in  order  to  model  the 
physical  effects  of  acoustic  energy  interacting  with  bone,  soft  tissue  and  hearing 
protection  devices.  The  original  method  is  a  pressure-based  formulation  of  the  acoustic 
wave  equation.  The  pressure-based  method  simulates  only  sound  wave  propagation 
through  inhomogeneous  fluid  media.  The  second  method  is  a  displacement-based 
approach  that  simulates  elastodynamic  wave  propagation  through  a  solid  medium  and/or 
a  fluid-solid  interface.  Both  algorithms  are  based  on  a  3D  finite  element  (FE)  matrix 
formulation  of  the  acoustic  wave  equation.  Novel  aspects  of  the  implementation  include 
the  fact  that  the  program  accepts  stacks  of  digital  images  to  form  the  3D  input  volume, 
the  program  uses  Message  Passing  Interface  (MPI)  directives  in  order  to  run  on  an 
intensively  parallel  cluster,  the  program  uses  a  uniform  grid  and  has  special  memory 
saving  techniques  making  it  possible  to  simulate  relatively  large  areas.  The  program  also 
uses  an  explicit  time-domain  algorithm,  such  that  no  matrix  inversion  is  required,  and 
multiple  frequencies  can  be  evaluated.  Also  a  Mur-Engquidst  absorbing  boundary 
condition  is  imposed  on  the  outer  boundaries  of  the  simulation  area  so  that  the  modeling 
of  radiation  and  scattering  effects  from  the  skull  are  more  accurate.  The  pressure-based 
algorithm  is  used  to  simulate  the  corresponding  experiments  on  live  human  subjects 
whereby  sound  is  transmitted  though  the  skull  via  head  mounted  vibrators. 

2.1.  Theoretical  aspects  (Pressure-based  method) 

The  acoustic  wave  equation  is  given  as 
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and  mathematically  characterizes  the  pressure,  p,  as  it  evolves  in  both  space,  r,  and  time. 
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t,  due  to  input  pressures  on  boundaries  and  interfaces.  In  this  equation  p(r)  is  a  spatially 
dependent  density,  c(r)  is  the  sound  speed  and  D(r)  is  an  attenuation  coefficient.  This  is 
an  inhomogeneous  wave  equation  with  no  input  driving  force  on  the  right  hand  side.  All 
inputs  are  represented  by  pressures  on  the  surface  of  one  of  the  modeled  structures.  The 
Finite  element  discretization  converts  this  spatially  continuous  wave  equation  into  one  in 
which  the  pressure  is  described  at  discrete  points  in  space  indicated  by  p,  in  the  following 
relation 
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The  simulated  area  is  sliced  into  smaller  individual  elements.  For  the  3D  model 
these  elements  are  cubes  with  8  nodes,  one  at  each  comer.  The  nodes  are  numbered  0 
through  N-l  and  {ft}  is  a  vector  of  length  N  representing  the  pressure  at  all  the  nodes  in 
the  mesh.  For  the  initial  program  attenuation  is  not  included  and  the  inertial  term  time 
derivative  is  approximated  with  a  centered  difference  equation  such  that 

a2Pl  Pr-2P;  +  P;-* 

at2  At2  '  * 

Insertion  of  this  formula  into  equation  2  and  solving  for  {p,}”0  yields 

{p.},"-[MnKFsM  +  2{Pl},-{Pi},*A,-[QF].  (4) 

In  this  formula  both  the  stiffness  matrix  [Krtifr]  and  the  mass  matrix  [MJ  arc  N  x  N 
matrices  which  relate  the  node  i  to  every  other  node  in  the  mesh  while  [QJ  is  a  row 
vector.  These  matrices  are  defined  by 

Kps.-jJ-JVNjVNidV  M-^-jNjN.dV  QF  -  -|-J'NJ^-dS  (5) 

in  which  N,  and  N4  are  linear  interpolation  functions  defined  at  the  nodes.  The  numerical 
algorithm  implements  this  equation.  Because  first  order  brick  elements  in  the  FE  mesh 
arc  used  each  row  has  a  non-zero  value  for  only  the  ij  entry  for  which  j  is  one  of  the  26 
nodes  surrounding  i  resulting  in  sparse  matrices  for  stiffness  and  mass.  Through 
diagonalization  of  the  mass  matrix,  [M],  matrix  inversion  is  avoided  and  a  purely  explicit 
time-stepping  scheme  is  achieved. 

Through  the  use  of  benchmarking  techniques,  it  has  been  determined  that  at  least 
40  nodes  per  wavelength  are  required  to  accurately  model  a  pressure  wave.  This  means 
for  an  input  pressure  signal  centered  around  4  kHz  propagating  through  air  with  a  sound 
speed  333  m/s  the  minimum  wavelength,  assuming  a  bandwidth  of  2  kHz,  is  333/5000  = 
6.6  cm.  Therefore  approximately  4003  =  64  million  nodes  would  be  required  in  order  to 
simulate  an  area  60  x  60  x  60  cm  [0.60/(0.066/40>=400].  In  order  to  not  store  a  matrix  of 
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27  x  64  million  only  unique  rows  of  the  stiffness  matrix  are  formed  and  stored.  The 
process  is  described  in  more  detail  in  Section  2.4. 

2.2.  Theoretical  aspects  (Displacement-based  approach) 

The  pressure-based  formulation  does  not  account  for  shear  phenomena  found  in 
solids  such  as  bone,  helmets  or  other  protective  gear.  In  order  to  account  for  the  solid 
shear  effects  the  elastic  pressure  wave  within  solids  must  be  modeled  along  with  the 
fluid/structure  interface.  The  two  main  methods  for  extending  the  finite-element  code  to 
include  solid  shear  effects  is  to  used  a  pressure-based  equation  in  the  fluid  and  a 
displacement-based  method  in  the  solid  or  to  use  a  pure  displacement-based  method 
throughout  the  both  fluid  and  solid.  In  the  mixed  prcssure/displacement  formulation  the 
loading  between  the  fluid  and  solid  must  be  satisfied  by  enforcing  continuity  of  pressure 
and  displacement  along  all  interfaces  between  the  two  types  of  media.  The  disadvantage 
of  this  method  is  that  the  all  nodes  along  fluid-solid  interfaces  must  be  identified  and 
stored  and  a  separate  interface  matrix  must  be  assembled  and  solved.  In  pure 
displacement-based  equations,  particle  displacement  is  solved  for  both  the  fluid  and  solid 
media  and  therefore  continuity  conditions  at  fluid-solid  interfaces  are  automatically 
satisfied.  The  disadvantages  of  this  method  are  that  1)  because  displacement  is  a  vector 
quantity  3  degrees  of  freedom  must  be  stored  for  both  the  fluid  and  solid  and  2)  the 
irrotation  of  displacement  in  the  fluid  V  x  u  -  0 ,  is  not  automatically  ensured.  For  this 
project  the  displacement-based  method  is  adopted  as  it  is  much  easier  to  code,  the 
additional  memory  requirements  arc  incremental  for  the  high  performance  computers  and 
there  are  documented  methods  for  controlling  the  rotation  problem. 

Displacement  wave  propagation  through  solids  is  modeled  via  the  clastic  wave 
equation  written  succinctly  as  (Auld,  1990;  Eq.  6.5) 


V,.(C:V,u)-p^  (6) 

in  which  u  is  a  vector  of  the  displacements  in  the  x,  y  and  z  directions,  C  is  a  6  x  6 
matrix  of  the  elastic  coefficients  and  V,  is  a  strain  operator  written  out  as 
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If  one  were  to  expand  the  elastic  wave  equation  to  show  individual  terms  the 
result  would  be  three  equations  with  multiple  stiffness  terms.  For  instance  the  simplest 
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solid  is  the  cubic  type  for  which  there  arc  only  three  independent  values  for  the  elastic 
stiffness  constants  namely  cu,  c12  and  c44.  Expanded  the  wave  equation  appears  as 
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More  complex  solids  such  as  hexagonal  crystals  will  have  similar  forms  of  the 
above  equations  but  will  have  more  cross  terms  for  the  greater  number  of  independent 
stiffness  coefficients.  If  the  cubic  is  also  isotropic  then  cl2  -  cu  -  2c 44  and  there  are  only 

two  wave  velocities:  longitudinal  given  as  V,  -  -Jc^Tp  and  shear  given  as  V,  -  ,/c44  /  p  . 
For  solids  more  complex  than  isotropic  the  wave  velocities  are  direction  dependent. 

As  with  the  pressure  equation  a  finite  element  approximation  of  the  elastic  wave 
equation  yields  a  set  of  linear  ODE’s  that  can  be  expressed  in  matrix  format  as 
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Thus  there  are  nine  distinct  stiffness  matrices  such  that 
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Kxx--/vN.,CVN,ldV  K^.-j-VN^CVN^dV  - -/VN.CVN.dV  (10) 

Isotropic  solids  are  typically  characterized  by  Young’s  modulus,  E,  and  the 
Poisson  ratio,  y  •  Stiffness  coefficients  can  be  determined  via  the  following  relationships: 
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The  displacement  wave  equation  for  the  fluid  appears  as 
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which,  expanded,  yields  the  three  equations: 


(12a) 

(12b) 

(12c) 


in  which  B  -  pFc2  is  the  bulk  modulus  for  the  fluid  and  is  analogous  to  the  elastic 
stiffness  constants  in  the  solid.  These  equations  are  compatible  with  the  displacement 
wave  equations  in  solids  if  B-clt-cl2  and  cw=0.  Numerical  artifacts  and 
instabilities  may  arise  if  cu  -  0 ;  thus,  this  quantity  should  be  small  but  greater  than  0. 


Pressure  can  be  derived  from  the  displacement  via  p  -  -BV  •  ii  which  in 
numerical  form  becomes 
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23.  Preprocessing  of  images  and  mesh  creation 

The  3D  FE  mesh  is  created  from  digital  images  (see  Section  4).  The  desired 
simulation  is  of  sound  coupled  to  the  human  head  via  an  attached  transducer  at  different 
input  locations.  Digital  images  of  a  head  can  be  obtained  from  either  a  CT  or  MR]  scan. 
Typical  scans  of  cadavers  can  be  very  noisy.  If  one  looks  at  the  individual  pixel  levels  of 
these  images  there  are  many  different  shades  of  color  and/or  gray  level  values.  The  FE 
algorithm  requires  distinct  colors  or  pixel  values  in  which  each  value  corresponds  to  a 
material  type.  Crisper,  cleaner  images  were  obtained  by  CT  scanning.  The  scan  yielded 
360  2D  slices.  The  slices  were  imported  into  a  medical  rendering  software  package 
known  as  Amira  and  soft  tissue  representing  brain  fluid  and  skin  were  added  to  recreate  a 
replica  of  a  human  head.  Transducers  attached  to  the  mastoid  and  forehead  were  also 
drawn  in.  The  transducers  arc  included  in  the  simulation  as  initial  results  of  an  input 
signal  simulated  without  the  transducer  shows  the  majority  of  the  acoustic  energy  not 
being  coupled  into  the  head  but  radiating  away  from  the  skull.  Each  medium  in  the 
image  is  given  a  distinct  color.  This  replica  is  in  terms  of  acoustic  properties  such  that 
the  entire  brain  appears  as  a  homogenous  fluid-gel  material.  Examples  of  representative 
slices  arc  shown  below  for  both  the  original  scan  and  the  enhanced  drawing.  The  black 
(royal  blue)  area  is  background  and  has  the  acoustic  properties  of  air  for  the  simulation. 
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The  slices  are  exported  from  Amira  as  .tiff  files.  The  .tiff  files  arc  imported  into  Matlab 
and  a  simple  script  is  used  to  convert  the  colors  to  numbers.  For  instance  0  is  assigned  to 
black,  1  to  green  etc.  The  numbers  for  each  slice  arc  written  out.  in  straight  ASCII,  to  a 
.txt  file.  The  ,txt  files  are  numbered  in  sequential  order.  Additional  .txt  files  can  be 
added  at  the  beginning  and  end  of  the  sequential  list  in  order  to  have  more  background 
material  above  and  below  the  skull. 

The  program  reads  in  the  numbers  in  the  .txt  files  and  each  slice  is  assumed  to 
have  a  finite  thickness  equal  to  the  element  length  is  the  z  direction.  Thus  the  pixel 
values  for  the  first  slice  become  a  voxel  value  and  each  of  these  voxels  correspond  to  a 
cube  element  with  acoustic  material  properties.  These  material  properties  and  clement 
dimensions  arc  used  to  form  the  stiffness  and  mass  matrices  of  the  second  section. 
Separate  simulations  arc  run  for  input  signals  with  different  center  frequencies.  As  the 
frequency  of  the  input  signal  is  decreased  and  the  wavelength  increases  the  images  arc 
downsampled  so  that  a  distance  of  at  least  two  wavelengths  between  the  input  force  and 
the  outer  boundary  is  maintained.  In  Figure  1  are  cross-sectional  slices  showing  the  car 
canals.  In  Figure  2  are  cross-sectional  slices  showing  the  location  of  the  mastoid 
transducers.  Figures  I  and  2  show  different  cross  sectional  slices  of  the  same  3D  scan. 
Figure  2  is  a  lower  (caudad)  slice  than  Figure  I.  Figure  I  shows  the  development  form 
the  original  scan  to  the  final  version.  Figure  2  shows  3  different  slices  for  corresponding 
reductions  in  frequencies.  The  color  maps  in  the  2  figures  are  different.  Intensity  plots  of 
the  propagation  of  the  signal  with  input  at  the  mastoid  are  shown  in  Figure  3  for  8  kHz.  4 
kHz  and  I  kHz  center  frequencies. 


Figure  I:  The  original  slice  of  a  CT  scan  of  a  dry  skull,  the  enhanced  image  with  soft 
tissues  and  earplugs  added  (34  cm  wide  by  50.3  cm  long)  and  the  downsampled  image  to 
simulate  a  4  kHz  signal  (68  cm  wide  by  100.6  cm  long). 


2.4.  Array  structure  and  memory  conservation  techniques 

The  program  consists  of  4  major  sections  namely  the  input  section,  the  elemental 
matrix  assembly  section,  the  global  matrix  assembly  section,  and  the  time-stepping 
section.  The  input  routine  reads  in  the  digital  images  and  an  input  file  which  includes 
additional  information  about  the  model  such  as  transducer  location,  transducer  frequency, 
length  of  run  and  material  parameters.  Within  the  input  section  there  are  two  notable 
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routines.  One  is  constraint)  which  processes  input  information  to  determine  node 
numbers  for  excitation  inputs,  dement  size,  node  numbers  for  recording  time  traces.  The 
other  major  subroutine  of  the  input  section  is  ‘pointer’  which  reads  in  the  sequential  list 
of  digital  images  and  assigns  a  ‘node  type’  to  each  node  based  on  the  elements 
surrounding  that  node. 


Figure  2:  Slices  showing  the  placement  of  the  mastoid  transducers.  Original  (a:  far  left)  is 
34  cm  by  50.3  cm,  downsampled  by  a  factor  of  2  (b;  middle)  is  68  cm  by  100.6  cm. 
downsampled  by  a  factor  of  4  (c;  far  right)  is  136  cm  by  201.2  cm.  All  have  339  by  501 
pixels.  Each  pixel  is  one  element. 


Figure  3:  Intensity  plots  showing  acoustic  pressure  emanating  from  a  right  mastoid 
transducer  for  8  kHz,  4  kHz  and  I  kHz. 


Every  element  has  a  node  at  each  of  its  8  corners  and  every  (non-boundary)  node 
is  surrounded  by  8  elements.  The  pressure  at  each  node,  i,  is  determined  by  multiplying 
the  27  non-zero  entries  in  the  j*  row  of  the  matrix  lM|  'lKwrfJ.  The  nodes  are  sequentially 
numbered  0  through  N-l  starting  with  the  lower  left  node  of  the  bottom  slice.  The  rows 
are  not  unique.  For  instance  all  the  black  or  background  nodes  will  have  identical  values 
for  the  27  entries.  The  same  is  true  for  the  other  color  nodes.  Interface  nodes  between 
colors  may  also  be  repeated.  There  are  8  elements  (or  pixels)  surrounding  each  node.  If 
an  interface  node  between  two  materials  has  the  same  pattern  of  elements  around  it  as 
another  interface  node  then  these  2  nodes  are  of  the  same  type.  There  is  a  global  integer 
array,  of  length  N,  known  as  NodeType  and  for  each  node,  i,  NodcType[i)  identifies  the 


type  of  node.  An  array  known  as  PixType  identifies  the  type  of  pixel  or  element.  Node  i 
is  surrounded  by  PixType[i-l).  PixType(il.  PixType[i+NXl-IJ,  PixTypc[i+NXl], 
PixType(i+NLAY-l],PixType[i+NLAY],  PixTypc[i+NLA  Y+NX 1-11  and 

PixTypc[i+NLAY+NXl]  or  elements  A,  B,  C,  D,  E,  F,  G  and  H  in  Figure  4.  Note  that 
NX  I  is  the  number  of  nodes  in  a  row  in  the  x  direction  and  NLAY  is  equal  to  the  number 
of  nodes  in  a  slice  or  NLAY=NX  1 XN Y 1 . 


Relationship  between  node  and  elements  for 
global  stiffness  matrix. 


Node  numbering  for  elemental 
stiffness  matrix 


Figure  4.  Node  Elements  in  the  Computational  Model 


The  27  entries  of  a  [K*jfr]  row  are  stored  in  a  3  x  3  x  3  array  known  as  the  global 
stiffness  array  UNKF[i][jJ[k][IJ  (or  UNKF_AB[i][j][k][l]  in  which  A  and  B  are  X,  Y,  Z 
for  the  displacement-based  code).  The  fourth  dimension  references  the  node  type.  For 
instance  all  the  background  nodes  will  have  NodcTypc[i]=4  and  therefore 
UNKFli J(jJlkJ[4J  will  store  the  stiffness  array  entries  for  node  type  4.  The  number  of 
different  node  types  depends  on  the  number  of  different  materials  plus  however  many 
unique  interface  nodes  exist.  This  value  is  stored  in  a  variable  NumUniqucNodes  and  is 
always  an  integer  several  orders  of  magnitude  less  than  64  million. 

The  global  stiffness  array(s)  is  assembled,  for  each  node  type,  in  the  routine 
GlobMatAsscmb.  The  array  is  assembled  from  the  elemental  stiffness  array,  EKF[i]|j|, 
which  is  formed  in  the  routine  ElemMatAssemb.  The  elemental  stiffness  array  is  an  8  x  8 
matrix  and  relates  each  local  node  in  the  element  to  every  other  local  node.  See  Figure  4 
for  local  node  numbering.  There  is  a  third  dimension  in  EKF  that  references  the  material 
type  of  that  element.  Note  that  unlike  the  global  stiffness  matrix  for  the  nodes,  EKF  will 
only  have  the  number  of  types  equal  to  the  number  of  materials  in  the  model,  usually  3  or 
4,  since  elements  are  not  split  across  interfaces.  The  GlobMatAsscmb  routine  forms  the 
stiffness  matrix  by  mapping  EKF  from  the  local  node  into  UNKF  for  the  global  node 
type.  For  instance  referring  to  Figure  4  UNKF[0][0][0]  is  equal  to  EKF[6J10]  of  element 
A  while  UNKFHH0J10]  is  equal  to  EKF[6]U]  of  A  plus  EKF[7][0]  of  B.  An  array 
known  as  IntFaceNodcType  maps  the  element  type  of  EKF  to  the  node  type  for  UNKF. 
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The  same  general  principles  are  used  to  assemble  the  mass  matrix  array  but  it  is  much 
easier  since  the  mass  matrix  is  diagonal  according  to  a  well  established  FE  mass  lumping 
technique. 

2.5.  Parallelization  using  Message  Passing  Interface  (MPI) 

The  current  trend  in  supercomputer  architecture  is  towards  distributed  memory 
machines  or  clusters  in  which  larger  autonomous  processors  with  individual  memory 
banks  are  linked  together  through  a  high  speed  communications  port.  These  machines 
are  easier,  hence  cheaper,  to  build. 

In  order  to  take  advantage  of  availability  and  speed  of  distributed  memory 
machines  the  algorithm  is  adapted  to  these  clusters  using  MPI  routines.  In  a  shared 
memory  architecture  there  are  many  small  processors  but  they  all  access  the  same  core 
memory  and  there  is  a  main  interface  node  that  controls  all  the  other  processors,  whereas 
in  a  distributed  memory  machine  each  processor  has  its  own  memory  bank  and  all  the 
processors  work  independently  of  each  other.  Therefore,  there  needs  to  be 
communication  between  the  processors.  Programming  in  MPI  on  clusters  requires  a 
slightly  different  mindset  because  the  program  will  be  launched  on  all  the  processors 
when  it  is  executed.  Thus  when  a  program  is  running  on  a  10  node  cluster,  there  will 
actually  be  10  versions  of  that  program  running  and  10  values  of  all  the  variables  defined 
in  that  program.  It  is  possible  to  limit  statements  and  variables  to  an  individual  processor 
as  there  is  a  variable  known  as  my_jank  which  is  equal  to  the  processor  or  node  number. 
Note  that  processors  in  clusters  are  also  referred  to  as  nodes  and  should  not  be  confused 
with  the  nodes  in  the  finite  element  mesh.  Thus  to  define  a  variable  N=100  on  node  5  a 
program  could  have  the  statement  if(my_rank=5)  N=100.  Without  the  ‘if  statement 
there  would  be  10  variables  equal  to  100  on  all  the  nodes. 

The  current  program  is  adapted  to  clusters  by  assigning  one  slice,  in  the  digital 
image  set,  per  processor.  Thus  in  the  current  example  there  are  400  slices  with  400  x  400 
nodes  each  so  the  MPI  code  would  be  launched  on  400  processors.  The  slice  number  will 
correspond  to  the  processor  number.  Within  the  TimeStcpAlg  routine  the  inner  loop  will 
step  through  400  x  400  =  160,000  nodes  rather  than  64  million  nodes.  When  it  executes 
400  processors  will  step  through  160,000  nodes  rather  than  one  master  processor  stepping 
through  64  million  nodes.  Naturally  each  processor  will  need  to  know  also  the  pressures 
for  the  slices  immediately  above  and  below  it  or  the  z+1  and  z-1  slice.  This 
communication  is  done  with  MPI_Send  and  MPI_Recv  commands.  At  every  time  step 
before  the  pressures  are  computed  a  set  of  MPI_Send  and  MPI_Recv  commands  are 
given  so  that  the  processor  for  every  slice  has  the  information  from  the  adjoining  slice. 
Arrays  UFFr  and  UFBk  are  length  160,000  and  store  the  pressure  values  for  the  z+1  and 
z-1  slices  respectively.  It  is  important  to  have  MPI_Barrier  commands  before  and  after 
each  time-step  so  that  one  processor  does  not  get  ahead  of  (or  behind)  the  others. 

In  addition  to  the  routines  updating  the  pressure,  the  Pointer  subroutine  also 
requires  MPI  commands.  The  my_rank  variable  is  used  so  that  each  processor  reads  in 
only  one  slice.  As  with  most  of  the  variables  the  NodeType  array  will  be  defined  on  all 
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the  processors  and  will  have  a  length  of  160,000.  In  order  for  the  routine  to  determine 
what  type  a  node  is  the  pixel  or  element  types  from  the  adjoining  slice  must  be  used. 
Therefore  within  the  pointer  subroutine  the  pixel  type  information  from  the  adjoining  z-1 
slice  is  sent  to  the  z  processor  and  the  z  processor  uses  this  information  to  determine  node 
type. 


There  is  no  message  passing  among  the  nodes  in  the  ElemMatAssemb  and 
GIobMatAsscmb  routines.  Each  processor  will  have  its  own  value  for  NumUniqueNodcs 
and  its  own  set  of  EKF  and  UNKF  matrices.  This  means  there  will  be  duplication  of 
these  matrices  since  many  of  the  node  types  on  one  processor  will  be  the  same  as  on 
another  This  tradeoff  in  memory  conservation  is  justified  by  the  ease  and  simplicity  of 
having  each  processor  use  its  own  set  of  matrices  when  updating  the  pressures  at  nodes  in 
its  slice. 

2.6.  Computational  model  results 

Simulations  of  experimental  tests,  in  which  transducer  inputs  are  located  at 
different  head  locations  and  behavioral  assessment  of  sound  is  recorded,  are  used  to 
assess  the  efficacy  of  the  numerical  model.  Acoustic  pressure  level  measurements  in  the 
ear  (auditory  canal  and  cochlea)  and  back  of  the  head  for  transducers  at  the  right  mastoid 
(Table  1)  and  forehead  (Table  2)  arc  tabulated.  The  flexibility  of  numerical  simulations 
allows  one  to  determine  pressure  levels  at  any  point  in  the  model  area.  For  these  results 
the  earplugs  were  turned  off  by  setting  their  material  parameters  equal  to  those  of  air. 
Measurement  points  are  in  the  air  pathways  (auditory  canal)  of  both  ear  canals  and  in  the 
bone  connected  to  the  inner  ear  (cochlea)  as  well  as  in  the  bone  at  the  back  of  the  skull. 
The  fact  that  resulting  pressure  is  generally  higher  in  the  air  pathway  (auditory  canal) 
indicates  that  not  all  the  input  signal  is  being  coupled  to  the  skull  and  some  energy  is 
being  conducted  via  the  background  air  as  is  also  indicated  by  Figure  3.  More  complete 
coupling  is  being  achieved  with  the  experimental  setup.  It  is  possible  to  determine  the 
error  due  to  air  conducted  energy  by  modeling  the  transducer  without  the  head  target  and 
subtracting  the  total  incident  field  from  the  total  field. 

2.7.  Three-dimensional  simulation  results  of  the  occlusion  effect 

The  occlusion  effect  was  modeled  numerically  and  in  three  dimensions  by  giving 
the  earplugs  in  Figure  14  (Section  4)  a  sound  speed  of  1878  m/s  and  a  density  of  1200 
kg/m\  Plots  of  the  pressure  measured  in  the  right  car  canals  for  the  right  mastoid 
mounted  transducer  centered  at  3  kHz  are  shown  in  Figure  5.  Note  that  there  is 
considerably  more  ringing  in  the  occluded  car  and  the  magnitude  of  this  ringing  is 
frequency  dependent.  Results  showing  the  peak  of  the  ringing  to  the  input  are  displayed 
in  Table  3.  The  pressure  is  higher  in  the  occluded  car  than  in  the  unoccluded  for  all 
frequencies.  For  the  occluded  car  the  pressure  is  higher  at  1  kHz  and  4  kHz  than  at  2  and 
3  kHz. 


Table  I :  Acoustic  pressure  measurements  as  a  function  of  frequency  in  the  ears  at  the 
specified  locations  for  transducer  attached  to  the  right  mastoid.  The  pressure  level  (in  dB) 
is  that  recorded  at  the  indicated  site  (auditory  canal  or  cochlea)  relative  to  the  pressure 


cvcl  just  under  the  transd 

uccr  in  the  skin. 

- - - 

4  kHz 

3  kHz 

2  kHz 

1.5  kHz 

1  kHz 

Right  auditory  canal 

3.1  dB 

-1.0  dB 

-5.7  dB 

-21  dB 

-26  dB 

Right  cochlea 

-4.5  dB 

-4.0  dB 

-3.4  dB 

-18.5  dB 

-19.8  dB 

Left  auditory  canal 

-52  dB 

-51.6  dB 

-55.6  dB 

-70.5  dB 

-73.7  dB 

Left  cochlea 

-56  dB 

-56  dB 

-53  dB 

-22.5  dB 

-32  dB 

Back  of  head 

-77  dB 

-71.6  dB 

-67  dB 

-75  dB 

-71  dB 

Table  2:  Acoustic  pressure  measurements  as  a  function  of  frequency  in  the  ears  at  the 
specified  locations  for  transducer  attached  to  the  forehead.  The  pressure  level  (in  dB)  is 
that  recorded  at  the  indicated  site  (auditory  canal  or  cochlea)  relative  to  the  pressure  level 
just  under  the  transducer  in  the  skin. 


4  kHz 

3  kHz 

2  kHz 

1.5  kHz 

1  kHz 

Right  auditory  canal 

-52  dB 

-37.4  dB 

-54  dB 

-66  dB 

-76  dB 

Right  cochlea 

-59  dB 

-41.6  dB 

-55  dB 

-80  dB 

-70  dB 

Left  auditory  canal 

-60  dB 

-35.6  dB 

-54  dB 

-76  dB 

-79  dB 

Left  cochlea 

-70  dB 

-48.8  dB 

-66  dB 

-60  dB 

-79  dB 

Back  of  head 

-98  dB 

-74.8  dB 

-80  dB 

-105  dB 

-98  dB 

lAmur*  9%  ten*  lor  uwoedudod  car  SkH/ 


Pt*«*uf*  »«lm  for  occluded  »*r  JkHr 


*  >  «»  i  ii  t  a  i  a 


Tun*  (mS*c) 


Figure  5.  Plots  of  pressure  versus  time  in  the  car  canal  and  ear  canal  bone  for  the 
unoccluded  (left)  and  occluded  (right)  ear.  Input  signal  is  blue,  pressure  in  canal  air 
passage  is  green  and  pressure  in  canal  bone  is  red. 
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Table  3.  Simulation  results,  in  dB,  of  pressure  in  the  ear  canal  air  passage-way  as  a 
function  of  frequency  for  the  occluded  and  unoccluded  ear. _ _ 


1  kHz 

2  kHz 

3  kHz 

4  kHz 

unoccluded 

-19dB 

-19dB 

-24db 

-22dB 

occluded 

-9.5dB 

-1  IdB 

-15dB 

-5dB 

3.  CQNCENTRIC.SEHEREJ1N1TE  ELEMENT-MODEL 

Like  many  numerical  methods,  it  is  possible  to  have  divergent  solutions.  It  is 
therefore  important  to  validate  the  model  with  geometries  that  have  analytical  solutions. 
One  such  geometry  is  that  of  two  concentric  fluid  spheres.  One  can  approximate  the 
human  head  as  two  concentric  spheres  where  the  outer  sphere  has  the  bulk  fluid 
properties  of  bone  (in  this  case  neglecting  the  presence  of  shear  waves)  and  the  inner 
sphere  has  the  properties  of  water  (McNew  et  al.,  2009). 

The  choice  to  use  concentric  spheres  provides  a  simple  starting  point  for 
verification,  especially  because  of  the  availability  of  a  closed  form  solution.  Verification 
has  been  focused  on  the  comparison  of  specific  spatial  points  of  FEM  fluid  and  solid 
against  the  fluid  concentric  sphere  closed  form  solution.  These  points  have  in  addition 
been  studied  in  an  experimental  setting  (recorded  using  hydrophones,  microphones  and 
accelerometers).  Work  is  under  way  to  verify  the  model  results  against  these 
experimental  results. 

The  computed  ANSYS  FEM  models  were  based  on  a  3D  tetrahedral  mesh  with 
spherical  boundaries.  This  mesh  is  a  requirement  by  ANSYS  to  provide  infinite 
boundaries  (open  air).  Interpolation  inside  of  the  tetrahedra  was  used  to  find  the  pressure 
values  for  specific  spatial  locations.  The  Visualization  Toolkit  (VTK)  open  source 
library  was  used  to  facilitate  this  task. 

The  ANSYS  model  was  first  exported  into  a  file  format  compatible  with  VTK.  It 
was  then  read  and  analyzed  using  the  built  in  reader  and  probe  filter  in  VTK.  Probe  filter 
finds  the  interpolated  scalar  quantities  at  specified  locations  or  at  the  points  of  a  new 
mesh.  Pressure  waves  were  exported  into  a  format  readable  by  Matlab,  to  allow  for  quick 
computation  and  renderings. 

The  same  incident  signal  was  used  both  in  the  closed  form  solution  and  in  the 
FEM  models.  This  signal  is  a  raised  cosine  at  center  frequency  of  500  Hz  (Figure  6). 

y(t)- ^1  -  cos(^-j  jc°s(cot).  (14) 
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Figure  6.  Incident  pressure  wave. 


The  time  signals  shown  correspond  to  the  points  marked  in  Figure  7.  Time 
pressure  signals  superimposed  in  Figures  8  through  10  show  pictorially  the  level  of 
agreement  between  the  FEM  models  and  the  closed  form  solution.  The  plots  show  good 
agreement  between  the  different  solutions  with  amplitude  differences  between  solid  and 
fluid  simulations,  and  as  the  wave  propagates  into  the  concentric  spheres,  a  period 
difference  begins  to  appear.  Figure  8  shows  the  pressure  wave  at  an  early  stage,  still 
traveling  in  air.  Figure  9  shows  the  wave  at  the  center  of  the  spheres.  The  period 
difference  observed  at  the  center  is  AT  =  0.1  ms  or  a  shift  in  center  frequency  of  Af  =  37 


Ha. 


Direction  of  Propagation 


Figure  7.  Points  at  which  pressure  was  computed  (in  green). 
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x*-37.  y»-16,  z*-2 2cm 


Figure  8.  Pressure  at  reference  microphone  placed  outside  of  the  sphere  nearest  to  the 
acoustic  source. 


x=0,  y«0,  z=Ocm 


xa-7.5,  y«0.  z=Ocm 


Figure  10.  Pressure  near  edge  of  inner  sphere. 


3.1.  Comments 

A  finite  element  model  was  formulated  for  simulating  the  acoustic  wave 
propagation  into  concentric  spheres.  This  model  was  verified  against  the  closed  form 
solution  of  the  plane  wave  propagation  into  concentric  fluid  spheres.  Good  agreement  is 
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seen  between  the  analytical  and  numerical  solutions.  The  main  artifact  observed  in  the 
numerical  solution  is  a  small  shift  in  pulse  length  and  center  frequency. 

The  cause  of  this  numerical  inconsistency  can  be  related  to  the  mesh  size  (element 
edge  length  per  wavelength).  Figures  8,  9,  and  10  support  the  above  statement  by 
showing  how  the  error  becomes  more  prominent  as  the  plane  wave  propagates  in  the 
domain.  The  current  computational  limitations  make  it  difficult  to  obtain  finer  meshes 
(shorter  element  edges)  for  the  same  domain  size  and  medium  acoustic  properties, 
however  examining  smaller  problems  both  in  the  analytical  solution  and  FEM  solution 
could  be  a  feasible  compromise. 


4.  COMPUTED J1QMQG  R  APJULt  CTIIMAGE  SEGMENTATION 

To  obtain  a  3D  segmented  volume  from  Computed  Tomography  (CT)  images  of  a 
cadaver,  the  software  package  Amira  was  used.  DICOM  (Digital  Imaging  and 
Communications  in  Medicine)  slices  were  obtained  and  imported  into  Amira. 
Segmentation  was  performed  by  means  of  data  thresholding  (by  brightness  level)  and 
also  by  manual  volume  selection.  Amira  allows  for  complex  manipulations  of  the 
volume  to  be  segmented,  such  as  the  translations,  rotations,  growth,  addition,  and 
removal  of  material.  These  functions  became  essential  at  eliminating  artifacts  from  the 
CT  images. 

The  resulting  volume  (Figure  1 1)  contains  only  information  about  the  materials  of 
interest  present  in  the  object  (artifacts  have  been  removed).  Moreover,  each  material  is 
labeled  with  a  unique  identifier  (in  this  case,  a  unique  color).  This  allows  for  quick 
manipulation  at  a  later  time. 


Figure  1 1 .  Segmented  volume. 


Once  the  necessary  manipulations  have  been  performed,  different  materials  are 
segmented  into  layers.  These  layers  contain  non-overlapping  material  (no  piece  of 
material  can  be  in  two  or  more  layers).  Figure  12  shows  orthogonal  CT  cross  sections 
with  overlaid  material  boundaries,  and  a  3-D  rendering  of  the  segmented  volume. 


Figure  12.  Segmentation  with  Amira.  Top  left:  Segmented  volume  (transparency  effect 
added.)  Other  three  panels:  Orthogonal  CT  cross  sections  with  color  lines  showing  the 
boundaries  of  each  layer. 


Each  of  these  layers  is  assigned  a  color  (Figure  13),  which  is  used  for  display,  and 
in  our  case  as  the  input  to  a  routine  that  assigns  acoustic  properties  to  the  skin,  bone,  and 
brain. 
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Images  such  as  Figures  12  and  13  are  the  input  to  modeling  software  that  is  able 
to  assign  acoustic  properties,  and  signal  sources,  to  any  pixel  of  an  image.  The  addition 
of  earplugs  and  bone  conductors  to  the  model  meant  “drawing"  them  into  the  segmented 
volume  in  Amira. 

Bone  conductors  were  modeled  to  be  cylindrical  in  shape.  They  are  composed  of 
an  active  region,  were  the  signal  will  be  provided,  and  a  body,  representing  the  casing  of 
the  mechanical  parts.  A  perfect  coupling  between  the  active  region  and  the  skin  of  the 
model  was  assumed. 

Figure  14  shows  a  segmented  slice  of  a  dry  skull  CT-scan  (with  grown  in  skin).  In 
this  model,  Amira  was  used  to  draw  in  ear  plugs  (blue)  and  bone  conductors  (purple  and 
teal).  A  brighter  color  was  used  to  represent  the  active  area  (that  generate  sound)  for  the 
each  bone  conductor.  Ear  plugs  were  modeled  to  be  cylindrical  as  well.  They  have  a 
perfect  fit  with  the  ear  canal. 
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5.  ACOUSTIC  RAY  TRACING 

Acoustic  ray  tracing  is  a  technique  to  visualize  the  direction  of  propagation  of  an 
acoustic  wave  as  it  travels  through  a  medium.  The  idea  of  applying  ray  tracing  to 
acoustic  problems  has  also  been  studied  as  a  component  of  this  grant  (McNew,  2008; 
Beilina,  2009).  Herewith  is  proposed,  a  new  ray  tracing  approach  based  on  the  notion  of 
isosurfaccs  is  presented.  It  is  first  developed  for  a  uniform  medium,  and  a  plane  wave  is 
used  as  the  input  signal.  The  second  development  is  for  an  oblique  incident  plane  wave 
on  a  planar  interface  between  two  media  with  different  sound  speeds.  The  computed 
angles  are  compared  to  the  predicted  angles  by  Snell’s  law.  Ray  tracing  is  also  applied  to 
the  concentric  sphere  case. 

5.1,  Wavefront  reconstruction 

The  first  step  to  ray  tracing  from  a  scalar  field  is  to  reconstruct  wavefronts.  Pierce 
(1989)  describes  wavefronts  as  surfaces  in  which  a  given  feature  of  a  signal  arrives  at  the 
same  time.  Alternatively,  in  the  frequency  domain,  wavefronts  arc  defined  as  the 
surfaces  of  equal  phase. 
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Finding  wavefronts:  Another  way  of  phrasing  Pierce’s  definition  is  to  say  that  wc 
want  to  find  the  time  delays  between  an  incident  signal  ptlK(t),  and  each  point  of  the 

volume  p,yi(t)  in  question,  where  p,yi(t)  is  the  time  pressure  signal  obtained  by 

picking  point  (x.y.z).  We  can  of  course  perform  this  operation  by  means  of  cross- 
correlation: 

P«  (0  *  P».»,«  (0  -  /I Pi*  (t  +  •  (15) 

Alternatively,  we  can  perform  the  same  operation  by  finding  the  cross-correlation 
in  the  frequency  domain,  and  then  applying  the  inverse  Fourier  transform: 

F-,{F(ptnc(t)},F{p,y,(t)}}.  (16) 

where  F  is  the  Fourier  transform  and  *  is  the  complex  conjugate.  The  later  expression 
makes  more  efficient  use  of  computational  resources  because  the  Fourier  transform  of  the 
incident  signal  needs  only  to  be  obtained  once. 

Once  the  correlated  signal  is  computed  (p^),  the  lag  at  the  maximum  point  of 
correlation  is  found  by  taking  the  Hilbert  transform  of  pcwT ,  and  then  computing  the 
maximum  value  of  its  magnitude: 

max(|Hilbcrt(p  J).  (17) 

Matlab’s  command  max  will  return  both  the  maximum  value,  and  the  lag.  The  returned 
lag  is  then  stored  in  a  3D  matrix  with  the  same  number  of  points  as  the  original  volume. 
Therefore  the  steps  arc: 

1 .  Obtain  a  time  dependent  pressure  scalar  field  in  3D  volume. 

2.  Obtain  an  incident  waveform. 

3.  Take  the  ID  time  correlation  for  each  point  of  the  volume  with  the  incident 
signal. 

4.  Find  the  lag  corresponding  to  the  arrival  of  the  incident  waveform  at  that  point 

5.  Store  lags  in  an  arrival  time  volume. 


A  resulting  arrival  time  volume  (Figure  15)  contains  the  time  for  each  point  when  the 
incident  signal  arrived. 
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Figure  15.  Arrival  times  volume  from  a  plane  wave  traveling  in  a  uniform  medium. 


Finding  wavefronts:  From  an  arrival  times  volume  we  can  define  wavefronts  as 
isosurfaces,  that  is  surfaces  of  equal  lag  (Figure  16).  This  is  done  in  Matlab  simply  by 
the  isosurface  command,  which  returns  a  mesh  of  faces  and  vertices  corresponding  to 
elements  on  that  surface. 


Figure  16.  Reconstructed  wavefronts  from  a  plane  wave  traveling  in  a  uniform  medium. 
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Ray  Tracing:  Using  the  previously  generated  wavefronts,  the  next  step  is  to  draw 
the  path  that  a  signal  will  take  while  traveling  in  the  medium.  The  resulting  path  is  a  ray. 
The  general  procedure  for  obtaining  rays  from  wavefronts  is  as  follows: 

1.  Obtain  current  t,  and  next  isosurface  t2. 

2.  Find  normals  at  each  vertex  of  t, . 

3.  Compute  the  intersection  of  each  normal  with  the  planes  containing  each  face  of 
isosurface  t2. 

4.  Find  the  face  that  contains  the  intersected  point. 

5.  Set  the  intersected  points  in  isosurface  t2  as  the  new  starting  point  of  the 
algorithm  (i.e„  t, ). 

6.  Repeat. 

Figure  17  shows  a  snapshot  of  the  algorithm  between  two  isosurfaccs. 


Figure  17.  Triangles  from  two  isosurfaces  (t,  and  t2)  are  depicted  as  rays  normal  to  t, 
(starting  on  the  red  dots)  intersect  inside  the  triangles  on  t2  (ending  on  the  green  dots). 
The  process  is  repeated  for  the  next  isosurface. 


The  collection  of  rays  found  using  the  algorithm  outlined  above  is  drawn  together 
in  a  volume  to  show  the  direction  of  wave  propagation.  Figure  18  shows  the  rays  found 
in  a  uniform  medium.  Notice  that  all  of  the  rays  form  a  straight  line  along  the  direction 
of  propagation.  This  is  both  because  the  medium  is  uniform  and  a  plane  wave  stimulus  is 
used. 
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y  (dimensionless) 


Figure  18.  Rays  from  a  plane  wave  traveling  in  a  uniform  medium. 


5.2.  Plane  wave  at  oblique  incidence  on  a  planar  interface 

Results  were  obtained  on  a  domain  that  was  split  in  two  equal  rectangles.  The  first 
half  had  sound  speed  c, ,  and  the  other  c2 .  Other  physical  properties  were  kept  constant 
and  arc  not  included  in  this  study.  A  plane  wave  with  oblique  incidence  (not 
perpendicular  to  the  interface)  was  introduced  from  one  end  (F'igurc  19). 


P2* 

\ 

c2  \ 

C. 

CD 

Normal 

P> 

Figure  19.  Cross-section  of  domain  showing  two  domains  (green  and  blue)  with  different 
sound  speeds,  the  incident  and  transmitted  plane  waves,  and  their  respective  angles  with 
respect  to  the  normal  of  the  interface. 

This  problem  is  governed  by  the  ratio  of  sound  speeds  at  the  interface,  that  is,  by 
Snell's  law: 
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(18) 


c:  Sin(e2) 
c,  Sin(o,)  ’ 


where  0,  and  02  are  the  incident  and  transmission  angles,  respectively.  A  pressure  field 
was  computed  by  implementing  the  following  equations: 

p,  _  eXMCw(».)‘kirs-(<>.)) 

p2  -  c *(k > )♦ 1  jrsln(#3 ))  ^  (2Q) 

where  k,  and  k2  are  the  wave  numbers,  and  x  and  y  are  the  Cartesian  spatial  coordinates. 
Given  an  incident  angle  0, ,  we  can  use  Equation  18  to  calculate  02 : 


(21) 


Pressure  fields  were  obtained  from  Equations  19  and  20  as  a  function  of  sound 
speeds  and  incident  angle.  Arrival  time  volumes  and  wavefronts  were  reconstructed  from 
these  pressure  fields  for  a  variety  of  sound  speed  ratios  and  angles  of  incidence. 

In  contrast  with  the  uniform  field  arrival  time  volume  (Figure  15)  and  wavefronts 
(Figure  16),  the  arrival  time  volumes  and  wavefronts  obtained  here  show  a  line  of 
diffraction,  caused  by  a  change  in  sound  speed.  Figure  20  shows  such  an  arrival  time 
volume.  Figure  21  shows  the  resulting  wavefronts,  and  Figure  22  shows  the  rays  traced 
for  a  domain  with  c,  =  300  m/s,  c2  =  330  m/s,  and  an  incident  angle  0,  -  45° .  Line  x  = 
25  in  these  plots  shows  the  interface  between  the  two  halves. 

5.3.  Plane  wave  at  oblique  incidence  on  a  planar  interface:  Results 

Ray  tracing  simulations  were  run  for  different  sound  speed  contrast  ratios  and 
angles  of  incidence.  The  incident  angles  were  chosen  to  be  less  than  the  critical  angle  for 
each  sound  speed  ratio.  The  critical  angle  is  given  by  Snell’s  law  when  the  transmission 
angle  02  is  90*: 


(22) 


Each  arrival  time  volume  was  a  50x50x10  3D  grid.  The  predicted  angle  and  computed 
(mean)  angle  were  recorded  for  each  case  and  are  presented  in  Table  4.  There  was  full 
agreement  when  the  angle  of  incidence  was  0°,  as  expected.  Small  differences  are 
observed  between  the  predicted  and  computed  values. 
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Figure  20.  Arrival  time  volume  for  the  plane  wave  with  oblique  incidence  on  a  planar 
interface  case. 


i 


Figure  2!.  Reconstructed  wavefronts  from  the  oblique  incidence  on  a  planar  interface 
case. 


-25* 


Qi - 4 - »— - 1 —  —a.  ...  — 1 —  -  i  ...  . >.  _ _ 4 

0  §  10  IS  20  »  JO  35  40  45  SO 
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Figure  22.  Rays  traced  from  the  oblique  incidence  on  a  planar  interface  case  (top  view). 
The  dashed  line  shows  the  interface  between  the  two  domains. 

Table  4.  Comparison  of  predicted  incident  and  transmission  angles  (0,  and  0,)  and  the 
computed  angles  ( 6,  and  02 )  by  ray  tracing  for  three  sound  speed  contrast  ratios.  The 
computed  values  shown  arc  the  mean  angles  found  across  ray  segments.  E,  and  E2  are 
the  relative  errors  (represented  as  percent)  between  the  predicted  and  computed  angles, 

E,-|e,-e,|/el. 


c,  (m/s) 

Cj  (m/s) 

9, 

0. 

02 

02 

Bi 

E? 

300 

310 

0° 

0° 

0° 

0° 

0% 

0% 

10° 

10.12° 

10.34° 

10.33° 

1.2% 

0.1% 

30° 

30.05° 

31.11° 

31.10° 

0.17% 

0.03% 

45° 

44.97° 

46.94° 

46.88° 

0.07% 

0.34% 

60° 

59.97° 

63.49° 

63.19° 

0.05% 

0.47% 

300 

320 

0° 

0° 

0° 

0° 

0% 

0% 

10° 

10.17° 

10.67° 

10.70° 

1.7% 

0.28% 

O 

o 

30.11° 

32.23° 

32.22° 

0.37% 

0.03% 

45° 

44.91° 

48.96° 

48.39° 

0.2% 

1.16% 

< 

60° 

59.81° 

67.48° 

63.19° 

0.32% 

5.93% 

t 

300 

330 

0° 

0° 

0° 

0° 

0% 

0% 

10° 

10.13° 

11.01° 

11.03° 

1.3% 

0.18% 

30° 

30.06° 

33.37° 

33.20° 

0.2% 

0.51% 

45° 

44.93° 

51.06° 

50.94° 

0.16% 

0.24% 

60° 

59.76° 

72.29° 

72.09° 

0.4% 

0.28% 
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The  results  in  Table  4  suggest  that  the  ray  tracing  algorithm  is  accurately 
predicting  angles  of  incidence  and  transmission.  The  error  levels  are  for  the  most  part 
less  than  0.5%. 

An  interesting  trend  is  that  for  each  of  the  0,-10°  runs,  the  error  was  close  to 
1.5%  for  the  incident  angle.  Another  seemingly  different  error  estimate  is  found  for  the 
transmission  angle  for  the  300/310  m/s  sound  speed  ratio  and  60°  degrees  incidence  run. 
In  that  case,  the  error  was  5.93%.  It  is  believed  that  the  starting  conditions,  at  the  ray 
tracing  level,  arc  responsible  for  these  differences.  Ray  tracing  starts  with  the  choice  of 
the  first  isosurfacc,  and  then  a  number  of  isosurfaces  follow.  This  choice  is  not  the  same 
for  each  of  the  runs,  hence  the  difference  in  errors  between  incident  angles.  In  addition  to 
this  observation,  the  algorithm  docs  not  use  all  of  the  isosurfaces  found.  Instead,  a  down 
sampling  operator  is  used  to  choose  isosurfaces  that  at  multiples  of  a  number  (i.e.  every 
two  isosurfaces,  every  four,  etc).  This  is  done  in  order  to  control  computation  time. 

5.4.  Wavefront  reconstruction  applied  to  a  fluid  sphere 

A  low  contrast  simulation  was  run  to  test  the  ray  tracing  method  on  a  fluid  sphere 
submerged  in  a  uniform  medium.  A  3D  pressure  field  was  obtained  from  the  analytical 
wave  equation  solution.  The  model  parameters  for  the  fluid  sphere  are  listed  in  Table  5. 


Table  5.  Physical  properties  of  concentric  fluid  spheres  simulation. 


Medium 

t-i 

IS 

Sound  Speed 
(m/s) 

Radius  (cm) 

Outside  of 
Sphere 

l 

300 

Sphere 

1 

400 

9 

Figures  23  and  24  show  a  top  view  and  a  perspective  view,  respectively,  of  the 
fluid  sphere  ray  traces.  Figure  23  shows  that  rays  start  with  a  uniform  distribution  at  the 
source  of  acoustic  pressure  (y  =  50).  As  the  wave  propagates  towards  the  sphere,  rays  in 
its  vicinity  focus  into  it,  increasing  the  concentration  of  rays  at  the  center  of  the  image. 
Once  the  rays  leave  the  sphere,  their  distribution  is  no  longer  uniform.  It  is  evident  aty  = 
0  that  there  is  a  pattern  in  ray  concentration  for  the  rays  that  intersected  the  sphere. 

Because  this  is  a  3D  model,  the  acoustic  source  is  a  plane  wave,  and  the  scatterer 
is  a  sphere,  we  would  expect  that  the  focusing  effect  would  happen  symmetrically  in 
every  direction.  Figure  24  is  a  perspective  view  of  the  same  simulation  showing  the 
symmetry  described  above. 
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Figure  23.  (Top  view)  Rays  obtained  from  the  analytical  solution  for  a  water  sphere  (18 
cm  diameter.  5  mm  thickness)  in  air.  Incident  signal  was  a  raised  cosine  pulse  at  500  Hz. 
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Figure  24.  (Perspective  view)  Rays  obtained  from  the  analytical  solution  for  a  water 
sphere  (18  cm  diameter.  5  mm  thickness)  in  air.  Incident  signal  was  a  raised  cosine  pulse 
at  500  Hz. 
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5.5.  Comments 


A  ray  tracing  method  was  developed  to  visualize  the  acoustic  wave  propagation 
pathways  from  a  scalar  pressure  field.  This  technique  is  based  on  the  notion  of  equal 
pressure  wavefronts,  or  isosurfaces,  obtained  from  an  arrival  time  volume. 

To  find  wavefronts  from  a  scalar  field,  time  correlation  is  performed  at  each  point 
of  the  volume  between  the  incident  signal  and  the  total  acoustic  pressure  at  that  point. 
The  maximum  lag  found  in  this  operation  becomes  the  arrival  time  value  for  that  point  in 
space.  Having  found  the  wavefronts,  rays  are  traced  by  finding  the  intersection  of 
normals  from  one  wavefront  onto  another  wavefront  repeatedly. 

Verification  of  the  ray  tracing  technique  involves  generating  pressure  fields  with 
known  wavefronts,  and  therefore  with  known  ray  directions.  A  pressure  field  for  a  plane 
wave  incident  onto  a  uniform  field  was  found  and  the  ray  tracing  technique  was  used 
successfully,  finding  the  predicted  wave  propagation  direction.  A  plane  wave  with 
oblique  incidence  on  a  planar  interface  between  two  media  with  different  sound  speeds 
was  also  simulated,  and  the  rays  computed  were  verified  against  those  predicted  by 
Snell’s  law,  finding  good  agreement  between  them. 


6.  HUMAN  SUBJECT  TESTING 

Human  subject  testing  was  conducted  to  acquire  and  to  apply  empirical  data  from 
human  listeners  stimulated  with  bone  conduction  signals  to  develop  and  verify  a 
functional  model  of  bone  conduction  propagation  pathways  in  the  human  skull. 
Specifically,  the  aims  of  the  studies  were  (1)  to  evaluate  evidence  of  occlusion  effects 
from  bone  conduction  hearing  as  a  function  of  stimulus  location  and  frequency  in  open 
versus  closed  ear  canals,  and  (2)  to  determine  the  relation,  if  any,  between  static  force 
applied  to  the  bone-conduction  oscillator  and  hearing  sensitivity  (thresholds)  as  a 
function  of  oscillator  location  and  stimulus  frequency  in  open  ear  canals. 

All  participants  were  between  18-30  years  of  age  with  normal  hearing,  no 
structural  abnormalities  of  the  head  or  neck,  no  history  of  ear  surgery,  otologic  or 
neurologic  disease,  or  head  trauma. 

6.1.  Evidence  of  occlusion  effects  in  bone  conduction  hearing 

Purpose:  Evidence  of  bone-conducted  sound  pathways  through  the  skull  have 
been  explored  with  growth  of  loudness  tests,  auditory  brainstem  response  and  otoacoustic 
emissions;  however  a  comprehensive  understanding  of  these  pathways  has  not  been 
achieved.  Less  is  known  about  the  usefulness  of  behavioral  thresholds,  sound  pressure 
levels  in  the  ear  canal,  and  the  auditory  steady-state  response  (ASSR).  as  a  function  of  ear 
canal  occlusion  in  understanding  propagation  pathways  of  bone-conducted  signals  for 
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different  oscillator  placements.  Data  for  bone-conduction  hearing  arc  critical  to  test 
numerical  simulation  results. 

Procedure:  Behavioral  thresholds  (determined  in  1-dB  steps)  and  RMS  SPL 
measures  for  bone  conduction  stimuli  at  40  dB  HL  were  acquired  with  car  canals  open 
and  closed  in  four  young  adults.  Ear  canal  status  (open  vs.  closed)  and  bone  oscillator 
placement  (forehead  vs.  right  mastoid)  were  evaluated  across  the  frequency  range  of  250 
to  4000  Hz.  In  addition,  ASSR  measures  were  obtained  across  repeated  test  sessions  in  a 
single  subject  for  bone-conduction  stimuli  presented  at  40  dB  HL  at  four  frequencies:  500 
Hz*  1000  Hz*  2000  Hz,  and  4000  Hz  at  two  placement  locations:  forehead  and  mastoid. 

Result  Summary:  Thresholds  were  less  sensitive  in  the  open  condition  compared 
to  the  occluded  condition  from  250  Hz  to  2000  Hz.  In  addition,  the  RMS  in  the  ear  canal 
was  higher  in  the  occluded  conditions  below  2000  Hz.  This  trend  in  hearing  sensitivity 
and  RMS  values  were  reversed  above  2000  Hz.  ASSR  results*  a  suprathreshold  measure* 
were  also  lower  (more  sensitive)  in  the  occluded  condition  across  frequencies.  The 
ASSR  results  were  acquired  in  a  single-subject  across  multiple  repeated  tests  sessions. 
Sample  results*  shown  in  Figures  25  and  26,  provide  evidence  of  the  occlusion  effect  for 
both  behavioral  and  physiologic  measures  as  a  function  of  oscillator  placement.  Given 
that  the  ASSR  was  consistent  with  the  behavioral  results  and  did  not  offer  additional 
information  for  modeling  of  bone-conduction  propagation  pathways  in  the  human  skull, 
further  ASSR  testing  was  not  indicated  for  the  current  project. 
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Conclusions:  The  data  from  all  three  measures  are  consistent  with  the  occlusion 
effect,  that  is,  occluding  the  car  canals  modifies  bone  conduction  propagation. 
Simulations  and  functional  models  of  bone  conduction  propagation  must  account  for  the 
occlusion  phenomenon.  The  results  obtained  from  this  portion  of  the  project  provide 
representative  values  that  could  be  used  for  further  development  and  verification  (Brault 
et  al.,  2007:  Wismer  ct  al,  2008). 

6.2.  Behavioral  bone  conduction  thresholds  as  a  function  of  static  force,  oscillator 
placement  location,  and  head  size 

Purpose:  The  purpose  of  the  study  was  to  determine  whether  bone  conduction 
hearing  thresholds  differ  for  oscillator  placement  location  if  static  contact  force  is  held 
constant.  Thresholds  for  bone  conduction  hearing  are  determined  with  a  bone  oscillator 
using  a  standard  headband  designed  to  apply  5.4N  of  static  force  to  the  oscillator  (ANSI 
S3.6-1996  R2004)  in  clinical  settings.  In  actual  practice,  however,  the  static  contact  force 
varies  over  a  range  of  several  Newtons,  as  the  headband  is  stretched  to  accommodate 
head  size  and  oscillator  placement  location  (e.g.,  mastoid  or  forehead)  or  if  the  headband 
is  modified  for  use  with  young  children  or  for  comfort.  Typically,  thresholds  of  bone 
conduction  at  the  forehead  are  less  sensitive  than  those  at  the  mastoid,  however  the 
contact  force  to  couple  the  oscillator  to  the  head  at  these  placement  locations  differs. 
Therefore  the  relation,  if  any,  between  placement  location,  contact  force  and  hearing 
threshold  sensitivity  is  not  clear. 
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Procedure:  A  total  of  24  normal  hearing  participants  (12  male,  12  female)  took 
part  in  the  experiment.  Participants  were  separated  into  two  groups  (small  and  large) 
based  upon  head  size.  There  were  12  participants  in  the  small  head  group  and  12 
participants  in  the  large  head  group.  Head  size  measures  were  made  from  nasion  to  inion 
using  calipers  and  tape  measure.  Head  circumference  measures  were  made  using  the 
CDC  guidelines  and  the  International  10-20  system.  Behavioral  bone  conduction 
thresholds  were  obtained  in  1-dB  steps  from  250  Hz  to  8000  Hz  in  1/6  octave  bands  at 
three  placement  locations:  right  mastoid,  left  mastoid  and  forehead.  Three  static  force 
levels  were  evaluated;  2N,  5N,  and  variable.  The  variable  force  levels  were  produced  by 
the  Standard  clinical  bone-conduction  headband  (ANSI  S3.6-1996  R2004).  For  each 
participant,  the  static  force  levels  applied  by  the  standard  headband  were  measured.  The 
2N  and  5N  force  levels  were  produced  by  customized  headbands  which  had  been 
fabricated  and  calibrated  (See  Section  7)  to  provide  a  known  static  force  applied  to  the 
oscillator  when  coupled  to  the  head. 

Results  Summary:  Data  from  the  small  head  group  showed  the  mean  threshold  for 
variable  force  was  1.48  dB  HL,  for  2N  was  3.14  dB  HL,  for  5N  was  3.72  dB  HL.  Mean 
thresholds  by  placement  were  5.87  dB  HL  at  forehead,  1.51  dB  HL  at  left  mastoid,  and 
0.96  dB  HL  at  right  mastoid.  Mean  thresholds  by  frequency  showed  that  2000  Hz  (8.81 
dB  HL)  was  significantly  higher  than  other  frequencies.  There  was  no  effect  of  gender 
and  no  interaction. 

Data  from  the  large  head  group  showed  the  mean  threshold  for  variable  force  was 
0.51  dB  HL,  for  2N  was  -0.45  dB  HL.  for  5N  was  6.72  dB  HL.  Mean  thresholds  by 
placement  were  3.28  dB  HL  at  right  mastoid,  2.69  dB  HL  at  forehead,  and  0.81  dB  HL  at 
left  mastoid.  Mean  thresholds  by  frequency  showed  that  2000  Hz  (6.50  dB  HL)  was 
significantly  higher  than  other  frequencies.  There  was  no  effect  of  gender;  however  there 
was  an  interaction  between  force  by  head  size  and  an  interaction  between  placements  by 
head  size. 

The  group  of  graphs,  Figure  27,  display  the  interaction  between  force  by  head  size 
and  placement  by  head  size.  A  significant  interaction  was  found  for  the  large  head  group, 
but  not  for  the  small  head  group. 

The  second  group  of  graphs.  Figure  28.  displays  bone  conduction  behavioral 
threshold  as  a  function  of  force  for  the  small  and  large  head  listeners.  (Please  note: 
About  3%  of  all  thresholds  collected  were  at  the  limits  of  the  audiometer  (-10  dB  HL) 
and  thus  were  replaced  with  median  values  for  the  group.)  The  static  force  levels  were 
greater  for  the  large  head  group  than  the  small  head  group.  Forehead  placement  yielded 
greater  force  than  mastoid  placement  for  both  groups. 
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Conclusions:  Mean  behavioral  bone  conduction  thresholds  obtained  with  the 
standard  headband,  that  produced  variable  static  force  levels,  were  not  related  to  the  static 
force  level  in  the  small  or  large  head  size  group.  Forehead  oscillator  placement  yielded 
thresholds  that  were  less  sensitive  than  mastoid  placement.  This  is  in  agreement  with 
previous  research  (Dirks  and  Malmquist,  1969).  For  the  calibrated  force  levels  of  2N  and 
5N  the  small  head  group  did  not  show  a  significant  difference  between  mean  thresholds 
obtained  and  again  forehead  placement  yielded  the  least  sensitive  thresholds. 
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Figure  28.  Behavioral  thresholds  by  static  force  across  four  frequencies  as  a  function  of 
oscillator  placement  with  the  standard  headband  in  groups  of  adults  with  small  vs.  large 
head  size.  Group  One  Graph  Symbol  Key:  1  kHz  =  diamond  1.5  kHz  =  triangle  2  kHz  = 
X  3  kHz  =  circle  4  kHz  =  square 


These  results  indicated  that  the  standard  clinical  headband  used  for  bone 
conduction  testing  is  capable  of  producing  a  range  of  static  force  levels  that  is  dependent 
on  oscillator  placement  location  and  the  listener’s  head  size.  For  the  small  head  group, 
oscillator  placement  location  seems  to  have  an  affect  on  bone  conduction  behavioral 
thresholds,  although  static  force  does  not. 
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The  large  head  group  showed  a  different  pattern  of  results  for  the  calibrated  force 
levels.  Mean  behavioral  thresholds  at  5N  were  less  sensitive  than  those  obtained  at  2N. 
Right  mastoid  thresholds  were  less  sensitive  than  forehead  thresholds.  In  order  to  further 
investigate  the  different  pattern  of  results  obtained  for  the  large  head  group,  a  subset  of 
large  headed  listeners  were  re-tested  in  a  follow-up  pilot  experiment  (See  Section  6.3). 

6.3.  Behavioral  bone  conduction  thresholds  as  a  function  of  static  force  and 
oscillator  placement  location  for  large  heads  (Pilot  study) 

Purpose:  A  subset  of  5  large  headed  listeners  were  re-tested  and  the  repeatability 
of  the  results  obtained  in  the  initial  experiment  (See  Section  6.2)  were  thus  documented. 
In  addition  to  the  2N  and  5N  levels  of  static  force,  three  additional  calibrated  force  levels 
(3.5N,  6N  and  8N)  were  also  tested.  This  was  done  in  order  to  examine  the  pattern  (i.e. 
linear,  curvilinear,  etc.),  if  any,  between  force  and  threshold  in  large  headed  listeners. 

Procedure:  Five  listeners  (1  female,  4  male)  from  the  large  head  group  who  had 
participated  in  the  initial  experiment  (See  Section  6.2)  were  tested.  Behavioral  bone 
conduction  thresholds  were  obtained  in  1  dB  steps  from  250  Hz  to  8000  Hz  in  1/6  octave 
bands  at  three  placement  locations:  right  mastoid,  left  mastoid  and  forehead.  Five  static 
force  levels  were  evaluated;  2N,  3.5N  5N,  6.5N  and  8N. 

Results  Summary:  Data  from  the  re-test  of  5  large  headed  listeners  showed  mean 
thresholds  of  1.01  dB  HL  for  2N;  -0.12  dB  HL  for  3.5  N;  0.84  dB  HL  for  5  N;  -0.20  dB 
HL  for  6.5  N;  and  -0.28  dB  HL  for  8N.  Mean  thresholds  for  each  placement  were  4.40 
dB  HL  at  the  forehead,  -2.18  dB  HL  at  the  left  mastoid,  and  -2.47  dB  HL  at  the  right 
mastoid.  Mean  threshold  for  frequency  showed  that  4000  Hz  (-4.01  dB  HL)  was 
significantly  lower  than  other  frequencies. 

Conclusions:  The  data  obtained  from  the  re-test  of  5  large  headed  listeners 
showed  no  significant  difference  between  the  mean  behavioral  thresholds  obtained  with 
the  5  calibrated  static  force  levels.  Forehead  placement  yielded  the  least  sensitive 
thresholds.  These  results  arc  not  in  agreement  with  those  obtained  in  the  initial 
experiment  (Sec  Section  6.2),  however  fewer  participants  were  tested.  Additional  re¬ 
testing  of  large  headed  listeners  is  warranted  in  order  to  more  fully  determine  if  static 
force  and  head  size  are  related  in  listeners  with  large  heads  (Brault  et  al„  2008). 

6.4.  Summary  of  human  subject  testing 

Experimental  evidence  from  behavioral  tests  of  hearing  sensitivity  (threshold 
measures)  and  physiologic  responses  (supra-thrcshold  Auditory  Steady-State  Response, 
ASSR)  were  evaluated  with  the  ears  open  and  closed  (occluded).  The  ASSR  was  chosen 
instead  of  the  ABR  because  the  spectra  of  these  stimuli  used  to  obtain  the  ASSR  arc  more 
frequency  specific  than  those  used  to  obtain  auditory  brain  stem  response  (ABR) 
measures.  Further,  in  general  the  results  from  ASSR  are  in  excellent  agreement  with 
behavior  thresholds  for  both  high  and  low  frequencies.  In  addition,  its  acquisition  time  is 
considerably  shorter  than  ABR,  further  enhancing  its  accuracy.  In  addition,  the  RMS 
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sound  pressure  levels  measured  in  the  ear  canals  (Real  Ear  Measures)  were  evaluated 
with  ear  canals  open  and  closed.  It  was  not  necessary  to  purchase  a  new  Real  Ear  Probe- 
Tube  system  for  this  project,  because  a  system  currently  used  in  the  clinic  became 
available.  Measurement  setting  and  protocols  were  modified,  in  consultation  with  the 
manufacture’s  research  and  development  team,  to  use  bone-conducted  stimuli.  The 
comprehensive  data  were  used  to  validate  the  computational  model  of  bone  conduction 
propagation  in  the  skull.  Numerical  and  experimental  results  showed  that  the  occlusion 
effect  is  more  prevalent  at  low  frequencies  and  decreases  as  frequency  increases.  The 
results  indicated  that  the  occlusion  effect  is  nearly  absent  at  2  kHz  and  at  higher 
frequencies. 

Bone  conduction  hearing  sensitivity  is  greater  for  oscillator  placement  at  the 
mastoid  region  compared  to  the  forehead,  consistent  with  previous  research  (Dirks  and 
Malmquist,  1969).  On  average,  for  groups  of  individuals  with  small  head  size,  bone- 
conduction  hearing  sensitivity  was  independent  of  the  static  force  applied  to  the 
oscillator.  Results  for  groups  of  individuals  with  large  head  size  were  more  variable  with 
greater  individual  differences  and  in  test-rctest  reliability.  Pilot  results  for  a  small 
number  of  large  headed  adults  were  inconsistent  upon  retest  and  should  be  investigated  in 
future  work. 


7.  HEAD  BAND  FORCE  CALIBRATION  FOR  BONE  CONDUCTION  STUDIES 

A  set  of  headbands  were  fabricated  and  calibrated  for  bone  conduction 
experiments.  For  these  experiments,  it  was  necessary  to  maintain  a  known  bone 
conductor  contact  force  across  subjects,  and  different  measuring  points.  The  headbands 
were  made  to  produce  known  forces  for  varying  stretch  distances  d  (Figure  29).  A 
system  of  calipers  was  devised  to  quickly  obtain  the  distance  d  from  a  subject,  and  match 
it  with  the  headband  that  would  most  closely  provide  the  desired  force. 
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Figure  29.  Headband  calibration  setup. 
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For  each  headband,  a  force  versus  distance  relationship  was  obtained.  The 
headband  was  positioned  on  the  border  of  a  table  with  one  end  attached  to  a  variable 
weight.  The  force  was  obtained  using  Newton’s  second  law  F  =  mg  where  g  is  the 
acceleration  due  to  gravity  (9.81  m/s2).  A  weight  of  204  g  is  therefore  equivalent  to  2N 
of  force.  The  stretch  distance  d  as  measured  and  recorded. 

A  total  of  five  headbands  were  calibrated:  four  custom,  and  one  standard  (used  in 
Section  6  studies).  Force  measurements  between  2N  and  5N  were  obtained  for  each 
custom  headband,  and  between  2N  and  ION  for  the  standard  headband  (Figure  30). 

The  results  show  that  with  the  custom  headbands,  forces  between  2N  and  5N  can 
be  easily  achieved  at  a  generous  range  of  distances.  The  standard  headband  force  to 
distance  relation  shows  that  it  always  exerts  more  force  than  the  custom  head  bands  for  a 
given  distance. 
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Figure  30.  Force  vs.  distance  plots  for  four  custom  head  bands  and  the  standard  head 
band. 


8.  OVERALL  CONCLUSIONS 

The  highly  interdisciplinary  program  has  addressed  successfully  the  challenging 
technical  goal  set  forth  to  develop  and  validate  an  acoustic  wave  propagation  model  using 
well  understood  and  documented  computational  techniques  that  track  an  air-borne 
incident  acoustic  wave  propagated  around,  into  and  in  the  human  head. 

A  finite  clement  model  was  formulated  for  simulating  the  acoustic  wave 
propagation  into  concentric  spheres.  This  model  was  verified  against  the  closed  form 
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solution  of  the  plane  wave  propagation  into  concentric  fluid  spheres. 

A  ray  tracing  method  was  developed  and  verified  to  visualize  the  acoustic  wave 
propagation  pathways  from  a  scalar  pressure  field. 

Finally,  the  occlusion  effect  provides  a  relatively  robust  phenomenon  to  compare 
the  computational  findings  (Table  3)  with  the  behavioral  and  functional  findings  (Figures 
22  and  23).  The  contour  (or  difference  score)  between  the  open  and  occluded  car  is 
shown  in  Figure  31  and  shows  a  similar  trend  between  the  computational  findings  and  the 
behavioral  and  functional  findings.  The  similarity  is  more  clearly  shown  with  the 
normalized  (at  I  kHz)  difference  score  (Figure  32). 

In  conclusion,  the  computational  approach,  utilizing  fluid  elements,  provides  very 
good  agreement  with  human  subject  findings. 
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Figure  3 1 .  Contour  (or  difference  score)  as  a  function  of  frequency. 
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Figure  32.  Normalized  contour  (or  normalized  difference  score)  as  a  function  of 
frequency. 
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